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Modeling H2 formation in the turbulent ISM: Solenoidal 
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We present results from high-resolution three-dimensional simulations of the tur- 
bulent interstellar medium that study the influence of the nature of the turbulence on 
the formation of molecular hydrogen. We have examined both solenoidal (divergence- 
free) and compressive (curl-free) turbulent driving, and show that compressive driving 
leads to faster H2 formation, owing to the higher peak densities produced in the gas. 
The difference in the H2 formation rate can be as much as an order of magnitude at 
early times, but declines at later times as the highest density regions become fully 
molecular and stop contributing to the total H2 formation rate. We have also used our 
results to test a simple prescription suggested by Gnedin et al. (2009) for modeling 
the influence of unresolved density fluctuations on the H2 formation rate in large-scale 
simulations of the ISM. We find that this approach works well when the H2 fraction 
is small, but breaks down once the highest density gas becomes fully molecular. 
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1 INTRODUCTION 

All observed Galactic star formation takes place within 
dense, massive clouds of molecular gas known as giant 
molecular clouds (GMCs). Understanding how these clouds 
form and evolve is therefore a crucial part of the study 
of star formation on galactic scales. In the past, molec- 
ular clouds have been seen as quasi-static objects that 
form stars slowly over a long lifetime where the dynami- 
cal evolution of a cloud and the chemical evolution of the 
gas within it were only loosely coupled and were modeled 
separately. However, observations provide velocity disper- 
sions documenting the existence of sup ersonic random mo- 
tions on scales larger than ~0.1 pc (e.g. | Go ldrcich & Kwanl 
1974 IZuckerman fc Evans!ll974 lLarsonlll98il; lMverslll983l: 



tions I Ballesteros-Paredes. Hartmann fc Vazauez-Semadenil 



Perault et alj 19861: Solomon et al 



19921: lOssenkopf fc Mac Low! 1200 



al. 

:>2; 



Rom an-Duval et al . 2011: 
motions have been associated with 
bulence in the interstellar medium 
appreciation that GMCs are highly 
and that their formation and 



19871: i Falgarone et all 
Hever fc Brunt! 120041: 



Schneider et all l201ll h These 
compressible tur- 
(ISM) leading to 
inhomogeneous 
evolution are dom- 



inated by the effects of supersonic turbulent mo- 
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19991: Mac Low fc Klessenl 12004 lElmegreen fc Scald 12004 
Scalo fc Elmegreenll2004 V The dynamical evolution of the 
clouds is rapid, with a timescale comparable to those of the 
most important chemical processes such as the conversion of 
atomic to molecular hydrogen or the freeze-out of molecules 
onto the surfaces of interstellar dust grains. In this picture, 
the dynamics and chemistry of the gas are strongly cou- 
pled, with one directly influencing the evolution of the other, 
meaning that they must be modeled together. 

The main chemical constituent of the molecular gas is 
molecular hydrogen, H2, with other molecules such as CO 
being present only in small amounts, so in practice the study 
of the formation of molecular gas is usually simply the study 
of the formation of H2. The molecule forms in the interstellar 
medium primarily on the surface of dust grains. Its forma- 
tion in the gas phase by radiative association is highly for- 
bidden due to the molecule's lack of a per manent dipole mo- 
ment and occurs at a negligibly slow rate ()Gould fc Salpeten 
1963), while the gas phase formation via intermediate molec- 
ular ions such as H~ or Hj~ is strongly suppressed by the in- 
terstellar radiation field (Glover 2003) and cannot produce 
molecular fractions much higher than xn 2 ~ 10~ 3 . 

Given the relatively slow rate at which H2 forms, it 
is natural to ask whether it is possible to produce large 
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amounts of H2 quickly enough for a model involvi n g rapid 
cloud formation to be viable. Idover fc Mac Lowl ( |2007bh 
have shown that dynamical processes such as supersonic tur- 
bulence have a great impact on the effective H2 formation 
rate. The presence of turbulence dramatically reduces the 
time required to form large quantities of H2. The density 
compressions created by supersonic turbulence allow H2 to 
form rapidly, with large molecular fractions being produced 
after only 1-2 Myr, consistent with the timescale required 
by rapid cloud formation models. It is found that much 
of the H2 is formed in high density gas and then trans- 
ported to lower densit ies by the action of the turbulence 
|Federrath et al.ll2008al ). a phenomenon that certainly has a 
significant impact on the chemistry of the ISM. 

One issue not ad dressed in the study by 
iGlover fc Mac Low! (|2007bh was the sensitivity of these 
results to the nature of the turbulent velocity field. Most 
of the work that has been done to date on the numerical 
modeling of molecular cloud turbulence has focussed on 
either purely solenoidal (i.e. divergence-free) turbulence, 
or weakly compressive turbulence where the solenoidal 
modes d ominate over the compressive (curl-free) modes 
(see e.g. iKlessen. Heitsch fc Mac Lowl l200d; iKlessen) l200ll ; 



lOstriker, Stone fc Gammidl200ll; iLemaster fc Stonejl2008h . 
The study by IGlover fc Mac Lowl |2007bT ) is no exception, 
as it used the same setup for generat ing weakly compressive 
turb ulence as i n ear lier work by iMac Low et al.l (Il998l ) 
and IMac Lowl (|l999l V Recently, however, Federrath and 
collaborators have perform ed a number of s tudies o f fully 
compressive turbu l ence jFederrath et al.l l2008bl . 120091 : 
ISchmidt et~aT]|2009l : IFederrath et alj|20ich . They show that 
compressive turbulence produces a significantly broader 
spread of densities than solenoidal turbulence with the 
standard deviation of the density probability distribution 
functions (PDFs) differing by a factor of 3 at the same 
rms Mach number and argue that while solenoidal forc- 
ing of turbulence is likely to occur in quiescent regions 
with low star formation rates like in the Polaris Flare 
and or Maddalena's Cloud, regions with a higher star 
formation activity are more compatib l e wit h comp r essive 
turbulence (see also, IFederrath et all l20ld : iBruntl l20ld : 
IPrice. Federrath fc BruntlbOllI ). 

The influence of the wide spread of densities produced 
by compressively driven turbulence on the rate at which 
molecular hydrogen forms in the ISM has not previously 
been investigated, but given the strong density dependence 
of the H2 formation rate, it is plausible that the effect could 
be large. To address this issue, we have carried out a nu- 
merical investigation of the rate at which H2 forms in inter- 
stellar gas dominated by compressive turbulence, and how 
this compares to the H2 formation rate in gas dominated by 
solenoidal turbulence. 

The outline of our paper is as follows. In §2 we describe 
our numerical method, paying particular attention to the 
treatment of chemistry and cooling, as well as the method 
used to generate and maintain turbulence in the gas. In §3 
we present our results for the H2 formation rate, and discuss 
the distributions of density, temperature and H2 abundance 
generated in the simulations. We also use our results to test 
the sub-grid sc ale model for H2 form ation in turbulent gas 
put forward bv lGnedin et alj (j2009h . We close with a sum- 
mary of our findings in section 4. 



2 SIMULATIONS 
2.1 Numerical method 

2.1.1 Chemistry and cooling 

Modeling the thermal evolution of the gas in a meaningful 
fashion and having a full chemical model of the ISM can 
easily require one to track several hundred different atomic 
and molecular species involved in several thousand different 
reactions, even if reactions on grain surfaces are neglected 
(see e.g. t he UMIST Database for Astrochemistry, as de- 
scribed in IWoodall et all 120071 ). This is impractical to in- 
clude in a 3D hydrodynamical code, since it would have 
an extreme impact on the code's performance. In order to 
run time-dependent chemical networks efficiently alongside 
the dynamical evolution of the system one needs to select a 
number of chemical species and mutual reactions such that 
the chemical network can be solved in a short enough time 
while still a dequately describing the ov erall evolution of the 
system (see IGlover fc Mac LowH"2007al lbri. For our purposes 
we need to be able to follow the formation and destruction 
of H2 with a reasonable degree of accuracy. 

We have modified the FLASH v2.5 adaptive m esh re- 
finement code (|Frvxell et al.ll2000l : ICalder et aLll2002T l to in- 
clude a detailed atomic/molecular cooling function and a 
simplified but accurate treatment of the most important 
hydrogen chemistry. FLASH is a massively parallel code, 
developed by the Center for Astrophysical Thermonuclear 
Flashes at the University of Chicago. It has support for a 
variety of different physical processes, including magneto- 
hydrodynamical (MHD) flows and self-gravity. FLASH uses 
the PARAMESH library to manage a block-structured adap- 
tive grid and the Message-Passing Interface (MPI) for par- 
allelization. 

Our modifications add a limited treatment of non- 
equi librium chemistry treated in an operator-split fash- 
ion ici over fc Mac Lowll2007aH bh. During each hydro step, 
the coupled set of chemical rate equations for the 
fluid are solved using the impl icit integrator DVODE 
l|Brown. Byrne fc HindmarslJll989T ). together with the por- 
tion of the internal energy equation dealing with compres- 
sional and radiative heating and cooling, under the assump- 
tion that the other hydrodynamical variables (e.g. den- 
sity) remain fixed. The advection of the gas energy den- 
sity is handled as in the unmodified FLASH code. Chemi- 
cal abundances are tracked using FLASH'S standard tracer 
field implementation, and consistent multifluid advection 
IjPlewa fc Mullerlll999r ) is used to reduce the advection er- 
rors. 

By default, the internal energy in FLASH is computed 
by subtracting the specific kinetic energy from the total spe- 
cific energy, using the equation 



e = E- 



jvr 
2 



(i) 



where e is the specific internal energy, E is the specific total 
energy and v is the velocity. In regions where the kinetic 
energy greatly dominates the total energy due to truncation 
error this approach can lead to unphysical (e.g. negative) 
internal energies, giving inaccurate values for pressures and 
temperatures. This problem can be avoided by evolving the 
internal energy separately, using the equation 
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Table 1. Reactions in our non-equilibrium chemical model. 



Table 2. Processes included in our thermal model. 



No. 


Reaction 




Reference 


1 


H + H + grain 


— ► H 2 + grain 


1 


2 


H 2 + H -> H 4 


H + H 


2 


3 


H 2 + H 2 -> H 


+ H + H 2 


3 


■1 


H 2 + e" -> H 


f H + e" 


4 


5 


H + c.r. -> H+ 


+ e" 


See jgjj 


6 


H 2 + c.r. -> H 


+ H 


See EH 


7 


H 2 + c.r. -> H 


+ H+ +e" 


See Q 


8 


H + e" -> H+ 


+ e- +e" 


5 


9 


H+ + e" -> H 


+ 7 


6 


10 


H + + e + grain — > H + grain 


7 



Process 



Reference 



References: 

iMac Low fc Shulll <1986ft 



Hollenbach fc McKej fl97: 



Il979h 
andv 



. _ , 3: iMartin. Keogh fc Mand 

4: iTrevisan fc Tennyson! (1200211. 5: lAbel et all Jl997h 
iFerland et al.l l|l992l l. 7: IWeingartner fc Draind l|200ll l 



dps 



2: 

(1998). 
6: 



-jg + V ■ [(pe + P)v] — v ■ VP = 0, 



(2) 



where p is the density and P is the gas pressure. The method 
used within the FLASH code is determined via the runtime 
parameter eint_switch. If the internal energy is smaller than 
eintswitch times the kinetic energy, then the total energy 
is recomputed using the internal energy from Eq.[2]and the 
velocities from the momentum equation. We have found that 
by setting eintswitch = 10~ 4 , we are able to avoid any 
problems due to truncation error. 

We treat the cooling coming from metals by assuming 
that the carbon, oxygen and silicon in the gas remain in 
the form of C + , O and Si + , respectively , as in the previ- 
ous studies of I Glover fc Mac Low] (2007a b). In practice, in 
the absence of photodissociating radiation (see below), we 
would expect carbon and silicon to rapidly recombine, and 
for the carbon to be converted to CO once the H 2 frac- 
tion becomes large. Ho wever, we know from previous work 
jGlover fc Clarkll2011al lbh that the behaviour of the gas is 
not particularly sensitive to whether the dominant coolant 
is C + or CO. Cooling from C + alone can reduce the gas 
temperature to values around 15-20 K, and although CO 
cooling enables the gas to reach even lower temperatures 
(T ~ 10 K), in realistic models of GMCs, the character- 
istic temperature of the fully molecular gas is generally in 
the range of 10 - 20 K dGl over fc Clarkll20llbh . As the H 2 
formation rate does not have a strong dependence on tem- 
perature, the approximate nature of our thermal treatment 
will have little influence on the H 2 formation rate in the 
gas. However, making this simplification allows us to mini- 
mize the computational requirements for our simulations by 
using a considerably simplified chemistry that follows only 
four species: free electrons, H + , H, and H 2 . We follow di- 
rectly the fractional abundances of molecular hydrogen xh 2 
and ionized hydrogen x H + (where these symbols denote the 
fraction of the available hydrogen found in these forms) by 
adding to the FLASH code an extra field variable for the 
mass density of each species. The abundances of the other 
two species - atomic hydrogen (ih) and electrons (a; c ) - are 
computed from the two conservation laws: conservation of 
charge 



C+ fine structure cooling 
O fine structure cooling 
Si+ fine structure cooling 
H 2 rovibrational lines 
Gas-grain energy transfer 
Recombination on grains 
Atomic resonance lines 
H collisional ionization 
H 2 collisional dissociation 
H 2 formation on dust grains 
Cosmic ray ionization 



Glover & Mac Low 
Glover et al. (2010) 
Glover fc Mac Low (2007 ,. ! 
Glover fc Abel (2008) 
Holle nbach fc McKee 
Wolfire et al. (2003) 
Sutherland & Dopita 
Abel et al. (1997) 
See Tabled] 
Hollenbac h fc McKee 
Goldsmith & Langcr 




and conservation of the number of hydrogen nuclei 



IH = Xn,tot — Zh+ — X K 2 



= £ H + + Xq+ + X Si + 



(3) 



(4) 

where xu,tot is the total abundance of hydrogen nuclei in all 
forms, and x c + and x sl + are the abundances of ionized car- 
bon and silicon, respectively, which remain fixed throughout 
the simulations. These species undergo the reactions listed 
in Table Q] The radiative and chemical heating and cooling 
of the gas is modeled with a cooling function that contains 
contributions from the processes listed in Table [5] 

We have also m o dified our treatment of the adiabatic 
index 7. iBolev et al.l (|2007f ) have recently pointed out that 
as the temperature of molecular gas increases, its specific 
heat capacity at constant volume, c v , changes due to the 
fact that first the rotational and then the vibrational energy 
levels of H 2 become populated and that therefore c v can- 
not be considered constant and independent of temperature 
as has been often assumed in previous numerical studies of 
star formation. For this reason, we use a set of lookup tables 
constructed with the assumption that the H 2 ortho-to-para 
ratio has its thermal equilibrium value. In these tables, the 
specific internal energy e is tabulated as a function of tem- 
perature T and fractional abundance of H 2 (xa 2 ), T is tabu- 
lated as a function of e and xu 2 , an d the adiabatic index 7 is 
tabulated as a function of e (or T) and xu 2 ■ To compute the 
required values for 7 or convert from e to T (or vice versa), 
we interpolate between the values stored in the tables. 

To test our modified version of the FLASH code, we 
performed static and turbulent simulations using both our 
new FLASH im plementation and our exis ting ZEUS-MP 
implementation (|Glover fc Mac Low! l2007al lbl) of the same 
physics, and verified that the codes produced comparable 
results. 



2.1.2 Turbulent driving and hydrodynamics 

We have applied our chemistry model to simulations of 
forced supersonic turbulence driven by fully solenoidal 
(divergence-free or rotational ) and fully compressive (curl- 
free or dilatational) forcing (Fe derrath et alj|2008bl . 120091 . 
2010), as two limiting cases to investigate the influence of 
the nature of the driving on the formation of H 2 . These 
simulations use the piecew ise parabolic method (PPM) 
ijColella fc Woodwordlll984h implementation of the FLASH 
code to integrate the equations of hydrodynamics on 3D pe- 
riodic uniform grids with 256 3 grid points. 
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As a control parameter in our simulations, we use the 
rms velocity of the turbulence. We use this in preference to 
the rms Mach number because the latter quantity depends 
on the sound speed of the gas, and in our non-isothermal 
simulations this is not constant, but varies in both space 
and time. To excite a turbulent flow with a specified rms 
turbulent velocity, we include a forcing term f in the gas 
momentum equation 



We adopt a rate £h 



10' 



s 1 for the cosmic ray ioniza- 



<9v , VP „ 

_ + (v .V)v = -— +f. 



(5) 



We model the random correlated stochastic forcing term 
f such that it varies smoothly in space and time using 
the Ornstein-Uhlenbeck (OU) process. The OU process is a 
well-defined stochastic process with a finite autocorrelation 
timescale T. ft describes the evolution of the forcing term 
/ in Fourier space (fc-space) with the stochastic differential 
equation: 



d/(M) = fo(k)V c (k)dW(t) - f(k,t) 



T 



(6) 



where W(t) is a Wiener process, a random process that adds 
a Gaussian random increment to the vector field given in 
the previous time step dt, followed by the projection tensor 
V^(k) in Fourier space. The projection operator reads 



(7) 



vU k ) = <^( fc ) + (i - cw*) = c** + C 1 

where 5y is the Kronecker symbol, and V^j = Sij — 
kikj/k 2 and V\j = kikj/k 2 are the fully solenoidal and the 
fully compressive projection operators, res pectively (see e.g. 
ISchmidt et ai]|2009l : iFederrath etaHl2010h . 

By changing the value of the parameter £, we can de- 
termine the power of the compressive modes with respect to 
the total forcing power. For £ = 1 in the projection opera- 
tor, we obtain a purely solenoidal force field, and with f = 0, 
we obtain a purely compressive force field. Any combination 
of solenoidal and compressive modes can be constructed by 
choosing £ £ [0,1]. 

The large-scale stochastic forci ng that we use, as the 
one c l osest to the observa tional data (lOssenkopf fe Mac Low! 
120021 : iBrunt et al.| [2009). models the kinetic energy input 
from large-scale turbulent fluctuations, breaking up into 
smaller structures. We thus drive the modes k — [1, 3] in 
units of 221, where L is the box size. The forcing amplitude 
A(k) has a parabolic dependence on k, such that most power 
is injected at |fe| = 2 and A(l) = A(3) = 0. 



2.2 Initial conditions 

Using the forcing module described above, and starting from 
zero velocities, we excite turbulent motions in a box with 
256 3 grid points and of side length L = 20 pc, filled with 
initially uniform atomic gas, using periodic boundary condi- 
tions. We perform purely hydrodynamical simulations, and 
neglect any complications introduced by magnetic fields or 
the effects of self-gravity. Th e abundances for carbo n, oxygen 
and silicon were taken from lSembach et al.l (|2000f l and are: 
x c + = 1.41 x 10" 4 , x = 3.16 x 10" 4 and x Si + = 1.5 x 10" 5 . 
We assume that the dust-to-gas ratio has the standard solar 
value, and fix the dust temperature at 10 K in every run. 



tion of atomic hydrogen (reaction 5 in Table [TJ . In the case 
of molecular hydrogen, we assume that all of the ions 
produced in the reaction 

H 2 + c.r. -> Ut + e~ (8) 



are destroyed by dissociative recombination, yielding two 
hydrogen atoms, and so adopt a rate Ch 2 ,6 = 2.22£h for 
reaction 6 that includes this contribution as well as that 
coming from direct dissociation of the H2. For reaction 7, 
we adopt the rate £h 2 ,7 = 0.037£h- In both cases, we as- 
sume that the ratio between the H2 destruction rates and 
th e ionization rate of at omic hydrogen is the same as given 
in IWoodall et all |2007l ). 

We perform two sets of simulations with different initial 
number densities: no = 30 cm -3 and no = 300 cm -3 . For 
each initial density, we perform simulations with rms tur- 
bulent velocities of 0.4 km s _1 , 2 km s _1 or 4 km s _1 , and 
examine both purely solenoidal and purely compressive forc- 
ing in each case, meaning that we perform a total of twelve 
simulations. We evolve each simulation for ten dynamical 
times T = L/2v Ilns . For the first two dynamical times, the 
chemistry module is switched off, and the turbulence is al- 
lowed to reach a statistically stead y state l|Federrath et al.l 
l200Sl l2010l : iPrice fc Federrath|2O10h . After that, we consider 
the chemical evolution and follow the gas for a further eight 
dynamical times. Note also that in our later discussion of the 
time evolution of the H2 fraction, we take the time at which 
we switch on the chemistry module to be t — 0, meaning 
that the simulations run from t — —2T until t = 8T. 

For simplicity, we set the ambient radiation field 
strength to zero in all of our simulations, thereby avoiding 
the necessity of modeling the penetration of Lyman- Werner 
band photons into the simulation volume, and allowing us 
to focus purely on the influence of the turbulent density en- 
hancements on the overall H2 formation rate. We note that 
the mean column density through our low no simulations 
is approximately 20 Mq pc -2 , which is more than sufficient 
to adeq uately shield the H2 in the gas against photodisso- 
ciation (|Krumholz et alj|2009h . provided that the incident 
radiation field is close to t he standard Galactic valu e. We 
have shown in other work l|Glover fc Mac Lowll201ll ) that 
H2 formation in clouds with surface densities of this value 
or higher is primarily limited by the time required to form 
the H2, rather than by the influence of UV photodissocia- 
tion. We therefore would not expect this omission to have 
a large impact on our results. At late times, we will tend 
to under-predict the amount of atomic hydrogen in the gas, 
and to over-predict the amount of H2, particularly in our low 
density r uns, but previous work su ggests that the effect will 
be small (|Glover fc Mac Lowll201 J ). We note, however, that 
this approximation will break down for clouds immersed in 
UV radiation fields that are significantly stronger than the 
standard Galactic value (Glover, in preparation). 



2.3 Numerical resolution 

iGlover fc Mac Low! (|2007bl ) and iMac Low fc Gloverl ||2010h 
examined the sensitivity of the H2 formation timescale in 
simulations of the turbulent ISM to the numerical resolution 
of the simulation, using numerical resolutions ranging from 
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0.4 0.6 
crossing time 

Figure 1. Time evolution of the mass-weighted H2 abundance in 
units of the turbulent crossing time for no = 300 cm -3 and t) rms 
= 4 km s _1 in simulations with numerical resolutions of 128 3 grid 
points (dotted), 256 3 grid points (dashed) and 512 3 grid points 
(solid). 



64 3 to 512 3 zones. They found that there was some depen- 
dence on the numerical resolution of the simulation at early 
times, owing to the ability of the higher resolution to better 
model the details o f the highest density structures formed by 
the t urbulence (see lFederrath et al.ll2010l ; |Price fc Federrathl 
120101 ). although it should be noted that in these simulations 
the turbulence was not driven to a statistical steady-state 
before the switch-on of the chemistry, which will tend to ex- 
acerbate any resolution dependence. These previous studies 
found that although there remain some signs of resolution- 
dependence at 256 3 zones, the difference between the 128 3 , 
256 3 and 512 3 results is very small. However, these resolution 
tests were performed only for the case of solenoidal turbu- 
lence. Therefore, to test the sensitivity of H2 formation to 
numerical resolution in the simulations with compressively 
driven turbulence, we have performed a resolution study for 
the run with w rms = 4 km s _1 and no = 300 cm -3 . This 
is the run in which the highest densities are produced, and 
so if this is well-resolved, then it is reasonable to assume 
that our lower density and lower u rms runs will also be well- 
resolved. In our resolution study, we performed simulations 
with resolutions of 128 3 , 256 3 , and 512 3 grid cells. 

In Figure [T] we show how the mass-weighted mean 
abundance of H2 (defined in section 3.1 below) evolves in 
runs with different resolution during the first crossing time. 
We see that there is almost no difference in the evolution 
of the H2 abundance in the three simulations, and con- 
clude that a numerical resolution of 256 3 grid cells should 
be enough to accurately model the growth of the H2 fraction 
in our simulations. 



3 RESULTS 

3.1 Time dependence of H2 abundance 

To quantify the rate at which H2 forms in our simulation 
we compute the mass-weighted mean molecular fraction, 
(zh 2 )m, given by 



o.s 



nf 0.6 
1 0.4 



S 0.2 



Solenoidal forcing 



V 


= 0.4 km s" 1 


D 


= 2.0 km s 1 


U 


s = 4.0 km s~' 




).01 0.1 1 

time [Myr] 



0.X 



— 

■If 0, 



S 0.2 




Compressive forcing 



1.01 



0.1 



1 

time [Myr] 



10 



(lH 2 )j 



£i,j,fePH 2 fej, k)AV(i,j, fc) 
M H 



(9) 



Figure 2. Evolution with time of the mass-weighted mean H2 
fraction (xjj 2 )m in runs with mean densities of 30 cm -3 (black) 
and 300 cm -3 (red). Three different values of the rms turbu- 
lent velocity v IIns are considered: 0.4 km s — 1 (dotted), 2 km s _1 
(dashed) and 4 km s — 1 (solid). The upper panel shows the re- 
sults for purely solenoidal forcing, while the lower panel shows 
the results for purely compressive forcing. 



where we sum over all grid cells, and where pH 2 (i, j, k) is the 
mass density of H2 in computational cell (i, j, fc), AV(i, j, fc) 
is the volume of the cell (i,j,k), Mn is the total mass of 
hydrogen present in the simulation. In Figure [2l we plot the 
evolution of (ih 2 )m as a function of time for both sets of 
runs, comparing different mean densities, rms velocities and 
types of driving. In Tabled we give the time in Myr required 
for the mass- weighted mean molecular fraction to reach 50% 
(t 50% ) and 90% (t 90% ). 

Looking at the evolution of H2 fraction with time in 
Figure [2] we see that the time required to convert a large 
fraction of the initial atomic hydrogen to molecular hydro- 
gen decreases as we increase the density or the strength 
of turbulent drivin g , in lin e with the previous findings of 
iGlover fc Mac Low! (|2007bl ). Comparing the two panels, we 
see that compressively-driven turbulence leads to more rapid 
formation of H2 than turbulence driven by solenoidal forc- 
ing. The difference is particularly pronounced at early times, 
and in runs with high rms velocities: for instance, t 50 % is 
roughly a factor of ten smaller in the compressive run with 
u lms =4 km s _1 and no = 300 cm -3 than in the correspond- 
ing solenoidal run. At later times, the differences between 
the compressive and solenoidal runs become much smaller, 
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Figure 3. As Figure [2] but showing the evolution of (xh 2 )m as a function of the turbulent crossing time T, rather than the absolute 
time. The left-hand panels show the evolution of (rrj^M from t = to t = 4T, while the right-hand panels zoom in on the period 
between t = and t = 0.5T. As before, we plot results for three different values of the rms turbulent velocity - 0.4 km s — 1 (dotted), 2 
km s — 1 (dashed) and 4 km s — 1 (solid) - and two different mean densities — 30 cm" " 3 (black) and 300 cm" 3 (red). 



with t 90 % varying by less than a factor of three even in the 
most turbulent runs. 

In Figure O we show the evolution of the mass- weighted 
mean H2 abundance as a function of the turbulent crossing 
time. Here we see that most of the dependence on the rms 
velocities vanishes when the time is measured in units of the 
crossing time. Regardless of the strength of the turbulence 
or the nature of the forcing, the molecular fraction reaches 
50% within only 0.1 - 0.2 crossing times in the high density 
model. For the low density case it takes approximately 0.5 
-1.0 crossing times to form the same amount of molecular 
gas, regardless of t) lms . 

Larger rms velocities yield more dense gas, resulting in 
a broader density PDF. On the other hand, they also lead to 
shorter turbulent crossing times, leaving less time for H2 to 
form. As shown in Figure[3] these two effects largely compen- 
sate for each other. In the solenoidal case, the latter effect 
dominates, and the H2 formation timescale, in units of the 
crossing time, decreases with decreasing u rma . In runs with 
compressive forcing, on the other hand, the increased width 
of the density PDF with increasing v Ilns is the dominant 
effect. 



3.2 Density and temperature distributions 

As Table [3] demonstrates, the H2 formation time does not 
scale linearly with changes in the density of the gas. We 
find that an increase in density by a factor of ten causes 
the gas to become 90% molecular only 5-8 times faster in 



Table 3. Time in Myr when the gas becomes 50% and 90% molec- 
ular in all our runs. 



Initial number density 


no = 30 


cm 


no = 


300 cm -3 


Solenoidal forcing 


*50% 


*90% 


t 50% 


*90% 


u rms = 0.4 km s 


17.94 


60.97 


1.91 


7.36 


D rms = 2.0 km s — 1 


4.79 


15.30 


0.64 


2.96 


Urms = 4.0 km s 


2.88 


9.67 


0.38 


1.83 


Compresive forcing 


*50% 


*90% 


*50% 


*90% 


Dims = 0.4 km s 


10.95 


42.73 


0.9 


6.74 


D ms = 2.0 km s — 1 


0.87 


6.74 


0.11 


1.44 


Urms = 4.0 km s 


0.36 


3.73 


0.036 


0.74 



the solenoidal case and 4-6 times faster in the compressive 
case for the same rms turbulent velocities. The reason we 
see less dependence than one might naively expect is clear if 
we look at how the density distribution varies as we change 
the mean density no- In Figure [3] we plot a volume- weighted 
number density PDF at t = 0.5 crossing times. As we de- 
crease the density, the entire PDF moves to low densities. 
Most of the H2 forms in dense gas, and so it is not surpris- 
ing that reducing the amount of dense gas available has a 
significant effect on xn 2 - However, the densest gas quickly 
becomes fully molecular and thereafter does not contribute 
to the total H2 formation rate (see Fig. [5} , reducing the ef- 
fect of density increase on the amount of formed H2. We 
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Figure 4. Volume-weighted density PDF for solenoidal (top) and 
compressive (bottom) forcing at time t = 0.5 crossing time in runs 
with v TulB = 2 km s _1 . The red solid line presents the PDF in 
the run with mean density 300 cm -3 , while the black dashed line 
shows the PDF in the run with mean density 30 cm -3 . 



therefore find a smaller difference between the H2 formation 
rates in the solenoidal and compressive runs than one might 
expect given the significant difference in the density PDF. 

In order to more quantitatively describe the H2 distri- 
bution, we examine how the H2 fraction varies with density. 
We compute xh 2 and n for each of the cells in the simula- 
tion volume and then bin the data by number density. We 
then compute the mean and standard deviation for xu 2 m 
each bin. The resulting values at t = 0.5 turbulent crossing 
times after the chemistry module is turned on are plotted 
in Figure [S] We clearly see a considerable scatter in the 
value of xn 2 at a given density. However, there is still an 
obvious underlying trend in the distribution of xu 2 with n, 
telling us that high density gas is mor e highly molecular, as 
expected (e.g. iHollenbach et alJll97ll ). At this point in the 
high density simulation the gas is almost fully molecular, 
whereas in the low density case xn 2 — 0.3 for solenoidal and 
xh 2 — 0.7 for compressive forcing (see Fig. [3|. Despite this, 
however, there are regions where the H2 fraction is already 
much higher, and we can see that gas with a number density 
n > 10 3 cm -3 is already almost entirely molecular in all of 
the simulations. 

We also examine how the gas temperature varies as a 
function of number density in our simulations. Just as with 
the H2 fraction above, we use the temperature output from 



our runs, bin it by number density n, and then compute the 
mean temperature and the standard deviation in the mean 
for each bin. We plot the resulting values again at t — 0.5 
turbulent crossing time in Figure [6] Strong shocks present in 
the turbulent simulations lead to high post-shock tempera- 
tures that can reach several thousand Kelvin. In low density 
gas, these shocks cause a significant scatter in the temper- 
atures. In high density gas, their effect is less pronounced, 
owing to the significantly shorter cooling time. In the case of 
compressive forcing, the gas is found to have a wider range 
of densities than the gas in the case of solenoidal forcing. As 
discussed before, this is a result of the stronger compressions 
produced by the turbulent forcing. 

A final notable feature in the temperature distributions 
is the fact that in the low density solenoidal run, the tem- 
perature of the gas at logn > 3.5 is clearly higher than in 
the other runs. This occurs because in this run, there is still 
a significant quantity of atomic hydrogen present at these 
densities (see Fig. allowing heating due to H2 formation 
to contribute significantly to the thermal balance of the gas. 
In the other runs, the atomic hydrogen fraction at these den- 
sities is very much smaller, and H2 formation heating does 
not play a significant role in determining the gas tempera- 
ture. 



3.3 Dependence on the density clumping factor 

As we are using periodic boundary conditions in our simula- 
tions, which prevent any of the H2 molecules that form from 
escaping from the simulation volume, it is relatively straight- 
forward to show that the evolution of the mass-weighted 
mean H2 abundance with time is described by the following 
equation 



d(zH 2 )i> 
dt 



= (2Rn 2 (T, T d )x H n - D H2 x H2 n)i 



(10) 



where Rn 2 {T, Tj) is the rate coefficient for H2 formation 
on dust grains (reaction 1), and Du 2 is a destruction term 
depending on both temperature and density that accounts 
for the loss of H2 in reactions 3, 4, 6 and 7 in Table [T] In 
practice, the impact of this destruction term is very small, 
unless ih <C xh 2 , and so to a good approximation 



d(:TH 2 )t 
dt 



(2i?H 2 (r,r d )x H n)M. 



(ii) 



As it stands, Eq.[TT]is not particularly useful, as in order to 
solve for the time dependence of (ih 2 )m, we need to know 
how Rh 2 , Xh, and n are correlated, and how this correlation 
evolves with time. However, we can convert Equation II II to 
a more useful form if we make a few further approximations. 
First, when the fractional ionization of the gas is small, as 
it is throughout our simulations, we have xu — 1 — xh 2 , and 
hence 



d(:TH 2 )t 
dt 



(2Rn 2 (T,T d )(l - xn 2 )n)u. 



(12) 



Second, in our simulations we keep the dust temperature 
fixed, and we know that most of the gas has a tempera- 
ture that lies within the fairly narrow range of 10 - 40 K 
(see Fig. [HJ. As the dependence of _Rh 2 (T, Ta) on T is weak 
when the temperature is low, we do not introduce a large 
error by treating the gas temperature (and hence _Rh 2 ) as if 
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Figure 5. Mean H2 fraction, plotted as a function of the num- 
ber density n of the gas at time t = 0.5 crossing time in runs 
with u rms = 2 km s _1 that use solenoidal (top) and compressive 
(bottom) forcing. The red solid line indicates the runs with mean 
density no = 300 cm -3 , and the black dashed line indicates the 
runs with mean density no = 30 cm -3 . To compute these values, 
we binned the data by number density and computed the mean 
value of xn 2 for each bin. The standard deviation in the value of 
xn 2 in each bin is indicated by the error bars. 



it were uncorrelated with the density, allowing us to write 
Equation [12] as 



d(xH 2 )i> 
dt 



2R H2 ({T) M ,T d )((l - x H2 )n)i 



(13) 



where (T)m is the mass-weighted mean temperature. 

To proceed further, it is necessary to make an addi- 
tional assumption regarding the correlation between the H2 
fraction and the density. Given the presence of the turbu- 
lence, it is appealing to assume that this turbulence per- 
fectly mixes the gas on a timescale much shorter than the 
chemical timescale. If we make this assumption, then we can 
treat xn 2 as being uncorrelated with density, allowing us to 
rewrite Equation 1 131 as 



d(zH 2 )i> 
dt 



= 2Jfe 2 ((r>M,r d )((l-XH 2 ))M(n)M (14) 
= 2 J R H2 «r) M ,r d )(l - (ih 2 )m)C„(ti) v , (15) 



where (n)v is the volume- weighted mean of n, defined as 



(n) v = y J nAV - 



(16) 




number density log n [cm ] 



100 



3 



10. 



Compressive forcing 



n = 300 cm" - 
n = 30 cm 




12 3 4 

number density log n [cm'] 



Figure 6. Mean gas temperature plotted as a function of the 
number density n at time t = 0.5 crossing time in runs with 
u Ims = 2 km s~ x using solenoidal (top) and compressive (bottom) 
forcing. The red solid line indicates the run with mean density of 
300 cm - 3 and the black dashed line indicates the run with mean 
density of 30 cm -3 . The data were binned in a similar fashion as 
for Figure [3] The standard deviation in the mean value in each 
bin is also indicated. 



This quantity is related to the mass-weighted mean of n by 



{n)u = I pndV, 



1.4m H f a,.. 
— rr — / n oV, 

1.4m H , 2> 
lAn H (n)vV^ )V ' 
C n (n)v, 



(17) 

(18) 

(19) 
(20) 



where C n = (n 2 )v / (n)y is the density clumping factor, and 
where we have used the fact that p = 1.4tohW, and hence 
that M = {p)vV = lArn-n(n}vV . 

Equation [15] demonstrates that if our assumption of 
rapid mixing of the H2 were true, then the evolution of the 
mass-weighted mean H2 fraction in a gas cloud would be 
related in a very simple fashion to the mean density of the 
cloud a nd its density clump ing factor. This fact has been 
used by iGnedin et alj (|2009l ) as the basis of a simple sub- 
grid scale model of H2 formation for cosmological simula- 
tions, or for other large-scale simulations without sufficient 
resolution to model the small-scale structure within molec- 
ular clouds. They write the formation rate of H2 in a similar 
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crossing time crossing time 

Figure 7. Evolution of the effective clumping factor C n e g (lines) and the true clumping factor C n (symbols) as a function of the 
turbulent crossing time T in runs with mean densities of 30 cm -3 (black) and 300 cm -3 (red). We plot results for three different values 
of the rms turbulent velocity: 0.4 km s — 1 (top), 2 km s — 1 (middle) and 4 km s — 1 (bottom). The left-hand panels show the results for 
purely solenoidal forcing, while the right-hand panels show the results for purely compressive forcing. 
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form to Eq. I15| and argue that C n ~ 3 -10 in typical tur- 
bulent clouds. Gnedin fc Kravtsovl (|201ll ) further developed 
this idea, and showed that this sub-grid model does a good 
job of reproducing the dependence of the average atomic 
and molecular gas surface densities on the total hydrogen 
surface density that is observed in nearby spiral galaxies 
|Wong fc Blitzj|2002l ). and the dependence of the mean H2 
fraction on the t otal hydrogen column density observed in 
our own Galaxy l|Gillmon et alj|2006l ; IWolfire et al.ll200Sft . 

However, the fact that we see a clear correlation be- 
tween xh 2 and n in our simulations (see Fig. [5} implies that 
the assumption of rapid mixing that we used to derive Equa- 
tion is incorrect. In reality, it takes roughly one-third of a 
turbulent crossing time to fully mix material from overdense 
clumps into their lower density sur roundings for solenoidal 
turbulence (|Federrath et al.ll2008al ). and potentially longer 
than this for compressive turbulence. Therefore, a prescrip- 
tion such as that in Eq. [15] will overestimate the H2 forma- 
tion rate. 

Our present simulations of solenoidal and compressive 
turbulence provide a useful test-bed for quantifying the ex- 
tent to which Eq. 1151 and by extension the Gnedin et al. 
sub-grid model, overestimates the H2 formation rate. To do 
this, we define an 'effective' density clumping factor 



C T l 



d(:EH 2 >M/di 



2i?H 2 ((T) M ,r d )(l - (xH 2 ) M )(n>v' 



(21) 



and compute how it evolves with time in each of our simula- 
tions, using our results for (a;H 2 )M and (T)m discussed ear- 
lier. We then compare this with the true density clumping 
factor C n computed at a number of different times during 
the simulations. The results of this comparison are plotted in 
Figure [7] (which shows the evolution between and 4 cross- 
ing times) and Figure [5] (which shows an expanded view of 
the first 0.5 crossing times). 

We see that at the very earliest times in the runs, there 
is a reasonable level of agreement between our inferred ef- 
fective clumping factor C n , e s and the measured clumping 
factor C n - Our computed values of C n ,ea are typically some 
20-40% larger than C n , but an error of this magnitude is 
plausibly explained by our use of the mass-weighted mean 
temperature in our calculation of Rn 2 '■ in reality, the dense 
gas, whose contribution initially dominates the H2 formation 
rate, will generally be colder than this mean temperature. 

However, this initial level of agreement between C n , c s 
and C n is very quickly lost in most of the runs. In all of 
the simulations, the true clumping factor C n remains ap- 
proximately constant, varying by at most a factor of two 
in the compressive case, and by much less than this in the 
solenoidal case. On the other hand, in most of the runs, 
C n ,cS decreases rapidly with time; only in the low density 
solenoidal model it does remain approximately constant dur- 
ing the lifetime of the simulation. The strong and almost 
immediate decrease of the effective clumping factor visible 
in Figures [7| and [8] is caused by the increase in the H2 abun- 
dance in the dense gas. As the dense regions that initially 
dominate the H2 formation rate become almost fully molecu- 
lar, their contribution decreases rapidly, causing a significant 
fall in the mean H2 formation rate within the simulation, and 
hence a significant decrease in C n , c ff • This effect is particu- 
larly pronounced in the compressively-forced runs, owing to 
their broad density PDFs. If we closely compare the results 



plotted in Figure[5]with the time evolut ion of the H2 f r action 
shown in Figure [3] we can see that the iGnedin et alj (|2009T l 
approach starts to break down when the gas is about 30% 
molecular. In the high density solenoidal runs, the H2 for- 
mation rate is almost immediately overestimated by a factor 
of 2, while in the compressive runs, the rate is overestimated 
by a factor of 4 in the low density case, and by a factor of 
10 in the high density case. 

It is clear from this analysis that in most cases there is 
no simple way to relate the mean number density of the gas 
and the current mass-weighted mean H2 abundance to the 
current H2 formation rate, given the strong time variation 
that we see in C„, c ff • This time variation is absent only when 
the characteristic H2 formation timescale is longer than a 
turbulent crossing time, as is the case in our low-density 
solenoidal runs, as only in this case is our assumption of 
rapid turbulent mix ing justified. One mu st therefore be care- 
ful when using the IGnedin et al] <|2009l ) sub-grid model to 
describe the H2 formation rate in numerical simulations. 



4 SUMMARY 

We have presented the results of a study of H2 formation in 
the turbulent ISM that examines the influence of the ampli- 
tude and mode of both solenoidal and compressive turbulent 
driving. We have performed high-resolution 3D hydrody- 
namic simulations using the massively parallel code FLASH, 
which we have modified to include a detailed treatment of 
atomic/molecular cooling and the most important hydro- 
gen chemistry. Even though the chemical network we use is 
significantly simplified compared to the most detailed mod- 
els available, it performs with acceptable accuracy for our 
purposes. We have performed simulations with numerical 
resolutions of 128 3 , 256 3 and 512 3 zones, and have demon- 
strated that our results are well-converged in our 256 3 runs. 
Our results also serve as a proof-of-concept application for 
our implementation of our non-equilibrium chemical model 
within the FLASH adaptive mesh refinement code. 

We find that with both compressively and solenoidally 
driven turbulence, molecular hydrogen forms faster in gas 
with a higher mean density, or an environment with 
stronger turbulence. Although initially (during the first 
million years), H2 formation is significantly faster with 
compressive turbulence than with solenoidal turbulence, 
at later times the differences become smaller, with the 
time taken to reach a molecular hydrogen fraction of 
90% varying by at most a factor of three between the 
compressive and solenoidal runs. In almost all of our 
simulations, the gas becomes highly molecular within a 
much shorter time than the 10-20 Myr that would plau- 
sibly be required to assemble the cloud from the diffuse 
ISM ( Ballcstcr os-Paredes. Hartmann fc Vazquez-Semadenil 



ll999l ; lElmegreenll2000l ;l Hartmann et alj|200lf ). 

We have also shown that when time is measured in the 
units of turbulent crossing time, the H2 formation timescale 
becomes much less dependent on the strength of the tur- 
bulence. Increasing the strength of the turbulence produces 
more dense gas and reduces the time taken to form H2. How- 
ever, it also reduces the turbulent crossing time of the gas. 
In the solenoidal case, the reduction in the turbulent cross- 
ing time is the dominant effect, and so H2 formation takes 
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Figure 8. As Figure[7] but showing an expanded view of the first 0.5 crossing times. As before, three different values of the rms turbulent 
velocity « rms are considered: 0.4 km s _1 (top), 2 km s _1 (middle) and 4 km s _1 (bottom) - and two different mean densities - 30 cm -3 
(black) and 300 cm~ 3 (red). The left-hand panels show the results for purely solenoidal forcing, while the right-hand panels show the 
results for purely compressive forcing. 
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longer (in units of the crossing time) as we increase u rms . On 
the other hand, in the compressive case, the broadening of 
the density PDF is the dominant effect, and increasing w rms 
leads to a moderate decrease in the H2 formation timescale 
measured in units of the crossing time. 

The differences we have found between the compressive 
and solenoidal runs can largely be understood by consider- 
ing the differences in the density PDFs in Figure [4] Com- 
pressive forcing produces a much wider spread of densities 
than solenoidal forcing, and since the H2 formation rate per 
unit volume scales almost linearly with density when xh 2 
is small, this allows the compressive runs to form H2 much 
more rapidly at early times. However, rapid H2 formation in 
the dense gas leads to its conversion to fully molecular form, 
at which point it no longer contributes to the total H2 for- 
mation rate. This phenomenon occurs in both the solenoidal 
and the compressive runs, but has a greater effect in the 
compressive runs owing to the faster initial H2 formation 
rate in these runs. 

Finally, we have also u s ed th e results of our study to 
show that the lGnedin et all |2009T ) prescription for correct- 
ing for the influence of unresolved density fluctuations on 
the H2 formation rate in large-scale Galactic o r cosmological 
simula tions must be used with caution. The iGnedin et al.l 
(2009) prescription assumes rapid gas mixing, when in re- 
ality it takes about one-third of a turbulent crossing time 
to mix the material from overdense clumps into the low 
density regions in the case of solenoidal forcing, and pos- 
sibly eve n longer in the case o f compressively-driven tur- 
bulence (|Federrath e t al. 20083). We have shown that the 
effective clumping factor calculated with the assumption of 
rapid mixing over-predicts the H2 formation rate. In the case 
of high density and strong compressive forcing, the H2 for- 
mation rate can be overestimated by more than an order of 
magnitude at all but the very earliest times. For applica- 
tions where one simply wants to determine which regions of 
the ISM become H^-dominated (i.e. more than 50% molec- 
ular) and how quickly this occurs, their approach remains 
reasonably accurate, since C n ,eB shows little variation while 
{xh 2 )m remains small. On the other hand, if one is inter- 
ested in the final, e quilibrium state of the gas (as in e.g. 
iKrumholz fe Gnedirj201(il ). then this approach may be prob- 
lematic, as it will systematically over-predict the H2 forma- 
tion rate in highly molecular regions, with the result that 
the H2 abundance will reach equilibrium too rapidly. 
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